Discrete Analog of the Burgers Equation 
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We propose the set of coupled ordinary differential equations dnj/dt = n 



- Uj as a discrete 

analog of the classic Burgers equation. We focus on traveling waves and triangular waves, and 
find that these special solutions of the discrete system capture major features of their continuous 
counterpart. In particular, the propagation velocity of a traveling wave and the shape of a triangular 
wave match the continuous behavior. However, there are some subtle differences. For traveling 
waves, the propagating front can be extremely sharp as it exhibits double exponential decay. For 
triangular waves, there is an unexpected logarithmic shift in the location of the front. We establish 
these results using asymptotic analysis, heuristic arguments, and direct numerical integration. 

PACS numbers: 05.45.-a, 02.60.-Lj, 02.50.-r 



o 

CD 



T3 

o 
o 



>: 

O' 

p; 

ov 

o: 

(N 



X 



The classic Burgers equation 

n t + (n 2 ) x = vn x 



(1) 



is the simplest partial differential equation which incor- 
porates both nonlinear advection and diffusive spreading 
[lHfj| . This ubiquitous equation emerges naturally in the 
presence of dissipation, and it is broadly used to model 
traffic flows [l| , transport processes @, 0] , surface growth 
|, and large scale formation of matter in the Universe 

The Burgers equation has two important properties. 
The first is continuity: equation (JT|) can be written in 
the form n t + J x = 0, hence assuring mass conservation 
in the absence of sources or sinks. If we view the quan- 
tity n as a density, then the total mass is a conserved 
quantity, J^° oo dxn(x,t) — const., as long as the density 
vanishes, n(x) — > in the limits x —¥ ±oo. The second 
property is asymmetry: due to the nonlinear advection 
term, equation ([T]) is not invariant with respect to the 
inversion transformation x — > —x. 

Our goal is to construct a discrete (in space) counter- 
part of the Burgers equation that maintains these two 
properties. Thus we discretize the spatial coordinate but 
keep the time variable continuous, n{x, t) — > rij(t), where 
j is integer. The differential equation 



dt 



f( n ]-i) - f( n ]) 



(2) 



represents a continuity equation on a one-dimensional 
lattice. Indeed, a finite total mass M = YlljL-oc n i rc ~ 
mains constant, M = const., if two conditions are met: 
(i) a vanishing density nj -> as j -> ±oo, and (ii) 
a vanishing function / at the origin, /(0) = 0. To re- 
produce the nonlinear advection term in ([1]), we take a 
purely quadratic and positive function f(n) = n 2 . With 
this choice, we arrive at the set of nonlinear difference- 
differential equations 



drij 
~~dt 



This system of equations meets the two criteria of mass 
conservation and asymmetry. 

Immediately, we can point out an important differ- 
ence between the discrete equation (J3j> and the contin- 
uous equation ([T]). Let us treat the spatial variable in §3§ 
as continuous, j — > x, and replace the difference with a 
second order Taylor expansion. The result of these two 
steps is the continuous equation 



nt + (n 2 ) x = (nn x ) 3 



(4) 



By construction, the nonlinear advection term is the 
same as in (JTJ) . However, the viscosity equals the den- 
sity, v = n, whereas in the original Burgers equation, 
the viscosity is constant. We restrict our attention to 
positive densities, 



n > 0, 



(5) 



'i-i 



(3) 



so that the solutions of ([3]) are stable (avoiding a negative 
diffusion instability). We note that transport coefficients 
often depend on density or temperature; in fluid dynam- 
ics [3|, [j; for example, transport coefficients vary as y/T 
for hard-sphere gases. 

Equation ^ describes the evolution of the probabil- 
ity density in a two-body analog of the standard Poisson 
process [121 ]. For example, we mention a homophilic net- 
work growth process [Tj|. In the canonical random net- 
work model, a pair of nodes are chosen at random and 
subsequently, the two are connected by a link. This ele- 
mentary step is repeated indefinitely, and in finite time, a 
percolating network emerges. As a model of homophilic 
networks where only similar entities interact, we con- 
sidered the situation where only nodes with exactly the 
same degree can be connected [Tj]. The degree distribu- 
tion rij(t), that is, the fraction of nodes of degree j at 
time t, obeys the rate equation Q. The initial condi- 
tion rij(0) = Sjfi represents a disconnected set of nodes. 
In this network context, the quantity rij(t) is a proba- 
bility density, and mass conservation guarantees proper 
normalization, Ylj=o n i = ^' Moreover, the condition 
([5]) reflects that probability distribution functions are by 
definition positive. 
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FIG. 1: The traveling wave. Shown is the function rij(t) 
versus j at four different times: t = 110, 120, 130, 140. The 
results are from numerical integration of the Eq. subject 
to the step initial condition ([6]). 



In this paper, we discuss the solutions of the discrete 
equation ([3]) in view of the well known solutions of the 
continuous equation ([1]). Using a combination of theoret- 
ical and numerical methods, we analyzed the solutions of 
the discrete equation for the following standard initial 
conditions [![: (i) a step function resulting in a traveling 
wave, (ii) a localized delta function leading to a triangu- 
lar wave, and (iii) a complementary step function with an 
ensuing rarefaction wave. For all of these cases, we find 
that the discrete analog faithfully captures the primary 
features of the continuous Burgers equation. However, 
we also find subtle and interesting departures from the 
classical solutions in the first two cases. Therefore, in 
the rest of this paper, we focus on traveling waves and 
triangular waves. 



Traveling Waves 



We first consider the step function initial condition 



nj (0) = 



j<0 
j>0. 



(6) 



Using the convenient Adams-Bashforth method [HI, [lj| , 
we numerically integrate the rate equation ([3]) and find 
that the solution approaches a traveling wave (figure [T]) 



nj(t) -> G{j - vt), 



(7) 



in the long-time limit. The function G(z) specifics the 
form of the traveling wave and v is the propagation ve- 
locity. Of course, the function G(z) has the limiting be- 
haviors lim z _ i ._ 00 G(z) = 1 and lim 2 _^oo G(z) = 0. 

To find the propagation velocity, we first note that 
according to equations ([3J and flBJ), the density does not 
evolve, rij(t) = 1 when j < 0. Next, we define the mass 




FIG. 2: The negative- z tail of the function G(z). Shown is the 
quantity 1 — rij versus j — t at various times. The theoretical 
prediction (1121) with f3 = 1.59362 is also shown for reference. 



in the positive half space, M + (t) = YljLi n j(.t)> an d n °t e 
that summation of (J3J gives 



dM+ 
dt 



(8) 



Since M+(0) = and no(t) = 1, the mass equals time, 
M + = t. This fact, along with equation ([7]) and the 
limiting behaviors of the function G(z), gives the velocity 



1. 



(9) 



Hence, the propagation velocity, which is dictate d by 
mass conservation, agrees with the continuous value 17 1 . 

To characterize the shape of the traveling wave, we 
substitute the form into the governing equation ([3|. 
The function G(z) satisfies the nonlinear and nonlocal 
differential equation 



G'{z) = G 2 {z)-G 2 {z-l), 



(10) 



where prime denotes derivative with respect to z. We 
now use asymptotic analysis to obtain the leading asymp- 
totic behaviors of G(z) in the limits z — ► ±oo. 

Since G(z) — ¥ 1 when z —¥ — oo, we substitute 
G(z) = 1 — 4>{z) into equation (jlOp . and obtain a linear 
yet nonlocal equation for the correction function (j), 



^(*)=2fe(z)-#z-l)]. 



(11) 



Therefore, the correction decays exponentially, 
4>(z) ~ e^ z . By substituting this behavior into (jlll) . we 
find the decay constant f3 = 1.59362 as the nontrivial 
root of the transcendental equation f3 = 2(1 — e~^). As 
shown in figure [2] the numerical results confirm that 



1-G(«)~ e 



(12) 



when z — > — oo. This exponential behavior agrees, at 
least qualitatively, with the corresponding behavior in 
the continuous case flal. 
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FIG. 3: The positive-z tail of the traveling wave function 
G(z). Shown is 2 3 n<j versus 2 j ~ t at three different times. 



We now invoke the limiting behavior g 
and find 



1 as z 



9{z) 



Z < Zq, 
Z > Zq. 



-00 



(16) 



This family of solutions is parameterized by zq, a quan- 
tity that depends on the details of the initial conditions. 
Hence, the traveling wave is exponential everywhere be- 
hind the front which is located at t + zq. In particular, 
the negative- z tail is analogous to (fT2")l . Remarkably, the 
leading front of the traveling wave is perfectly sharp [20l - 
[22| as the function g vanishes beyond the front location! 

We view the double exponential decay (jT4j) as a dis- 
crete analog of a perfectly sharp front. Moreover, the 
effective viscosity v = n and the vanishing density ahead 
of the propagating front are together responsible for the 
extremely sharp front (l7j . 



In the complimentary z — > 00 limit, we expect that 
G(z) — > 0. Now, the positive term in (fTU)) is negligible, 
and the behavior is governed by the nonlinear and non- 
local equation G'{z) = -G 2 (z- 1). We use the WKB 
transformation [19( G(z) ~ VP exp(— ip) and arrive at the 
linear relation 



ij>{z) = 2ip(z- 1). 



(13) 



We thus obtain the exponential solution ip(z) — 7 2 Z and 
further, the prefactor ^(z) = (4 In 2)ip(z) . The constant 
7 cannot be determined in the realm of asymptotic anal- 
ysis; numerically, we obtain 7 = 0.818. Remarkably, the 
leading tail of the traveling wave is extremely sharp, as it 
follows the unusual double exponential decay (figure [3]) 



G(z) 



-7 2" 



(14) 



when z — > 00. Hence, the front of the propagating wave 
does not extend beyond a few lattice sites (see figure 
[T]). In the continuous case the traveling wave is specified 
by G con t(z) = i[l — tanh(z/2^)] and hence, the shape is 
symmetric as G cont (z) + G cont (— z) = 1. Therefore, the 
behavior in the discrete case differs in two ways. First, 
the shape of the traveling wave is asymmetric. Second, 
the leading tail of the wave follows a double-exponential 
decay, in contrast with the simple exponential decay in 
the continuous case. 

To gain further insight into the extremely sharp front 
(fT4| . we also study the continuous equation (|4]). We 
assume that the solution approaches a traveling wave, 
n(x,t) — > g(x — vt), with velocity v — 1 (again, the ve- 
locity is set by the continuity condition). From Eq. (QJ, 
the function g(z) obeys the ordinary differential equation 



g> - 2gg' + (gg')' = 0. 



(15) 



Integration of this equation is immediate, and using the 
limiting behavior g — ¥ as z — ¥ 00, we have g' = g — 1. 



Triangular Waves 

Under the Burgers equation {TJ), an initial condition 
with compact support necessarily evolves into a trian- 
gular wave. Moreover, the shape of the triangular wave 
is universal. Without loss of generality, we consider the 
localized initial condition 



nj(0) = 5 jt 



(17) 



As discussed above, according to the discrete equation 
([3]), the total mass is conserved. Moreover, there is no 
evolution in the negative half-space, rij(t) = for j < 0. 
Therefore, the mass in the positive half space equals 
unity, M + {t) =J2j>i = L 

Numerically, we integrate the discrete equation ([3]) 
starting with the initial condition (|17[) and find that the 
density adheres to the scaling form (figure [4j 



n 3 {t) 



—F ( — 



Further, the scaling function is triangular (figure 



(18) 



F(x) = 



x < 0; 

< x < 2; 



(19) 







x > 2. 



Therefore, the discrete equation ([3]) captures the basic 
features of the continuous equation (p}. 

The scaling behavior (fT8"l) -(fT9 j) follows from the in- 
viscid Burgers equation dn/dt + 2ndn/dj — 0. Indeed, 
mass conservation, in conjunction with the inviscid Burg- 
ers equation, dictates the scaling form (fT5| and the nor- 
malization J^°dxF(x) = 1. By substituting (fT5|) into 
(U} and ignoring the diffusion term, we see that function 
F obeys 



d 
dx 



F F 



0. 



(20) 




There are two solutions, F — x/2 and F = 0, and conse- 
quently, the scaling function is piecewise-linear as in (|19|) . 
Finally, the extent of the nonzero region, < x < 2, is 
set by the normalization condition. 

We now analyze the solution of the discrete equation 
([3]) in the asymptotic regime. The quantity n\ satisfies 



the closed equation dni/dt = 
condition ni(0) = 1, we have 



1 + t 



and using the initial 



(21) 



This density is inversely proportional to time, ni ~t~ l 
in the limit t — > oo. The structure of the rate equations 
([3|) suggests that in general, the densities rij are inversely 
proportional to time, 



A jt -\ 



(22) 



when t — > oo. Substitution of (j2"2")l into Q shows that 
this asymptotic behavior is compatible with the rate 
equation, and furthermore, yields a quadratic recursion 
relation for the prefactors, 



Aj - A, 



A 2 



(23) 



Starting with Ai = 1, we can iteratively obtain the pref- 
actors from the explicit formula 



1 



leading to A 2 



_ 1 + V5 
2 



1 



AA)_ 



(24) 



A, = 



l + \/7+2v^5 



etc. I 25 



We use a continuum approach to calculate the large- j 
asymptotic behavior of Aj . To obtain the leading asymp- 
totic behavior, we convert equation (|23|) into the differen- 
tial equation dA/dj = 1/2. Hence, to leading order, the 
coefficients are linear, Aj ~ j/2, and from the asymptotic 



behavior (|22p one has rij ~ j/(2t). We thus recover the 
continuous behavior indicated by (fT5)l and F(x) = x/2. 
There is, however, a correction to the leading asymp- 



totic behavior. We substitute Ai 



3/2 



into 123 



and find that the correction Uj obeys the recursion 



(i-i)K' 



x.j_ x 



(25) 



and ui = 1/2. When j is large, the difference between 
the two quadratic terms is negligible, and the continuous 
approach gives du/dj = l/(4j). As a result, the prefac- 
tors include a logarithmic correction 



An 



Inj 



C. 



(26) 



The constant C — 0.32324 is computed by numerical 
iteration of flU) • 

The logarithmic correction affects the location of the 
propagating wave. If we denote the front location by j*, 
then mass conservation gives Yjj*=i tij — 1. Substitution 
of (|22|) and (|26[) into this sum yields 



J* 



2yft- 



(27) 



by straightforward asymptotic analysis. The leading 
term follows from the continuous behaviors (fT8| and the 
extent of the triangular wave specified by However, 
the logarithmic correction, which follows from the coef- 
ficient (|26p constitutes a departure from the continuous 
behavior. We conclude that the triangular wave has the 
following form 



/(*) 





Ajt- 1 



j<0, 

< j < j* 

j > j*- 



(28) 



We already established the front location and the 
maximal height of the wave n* ~ t~ x / 2 which follows 



5 



from (|T8|) and (jT9|l . Thus, we anticipate the scaling be- 
havior 

-i=GCj-j.), (29) 

in the finite neighborhood j — ~ 0(1)- Numerical in- 
tegration confirms this scaling behavior (figure [5]). The 
function G{z) specifies the shape of the front. By sub- 
stituting (|29|) and (|27l) into the rate equation ([3]), we 
see that the scaling function G(z) satisfies (fT0|) . Hence, 
triangular waves and the traveling waves have identical 
fronts! This behavior is remarkable in view of the two 
very different scaling forms (|T8]) and ([29]) . In particu- 
lar, the leading front is extremely sharp and follows the 
double exponential decay (TT4|) . 

We also studied rarefaction waves by considering the 
complementary step function 

mp)-{; : r.:; : (-) 

Again, there is no evolution in the negative half space. 
According to the continuous equation (TT]), the density 
obeys the scaling behavior rij(t) —> F(j/t) with the very 
same F(x) given by equation (|19j) . Of course, there is a 
significant difference with the triangular wave above as 
the extent of the wave is now linear in time. In addition, 
there are diffusive boundary layers, for example when 
j — 2i ~ O^ 1 / 2 ). Numerical integration of the rate equa- 
tion (j3|) shows that the scaling behavior and the diffusive 
boundary layer both agree with the continuous case. 

In contrast with the traveling and triangular waves 
above, we can replace the difference equation ([3]) with 
the differential equation (fTJ) even inside the front region 
because the size of the boundary layer grows diffusively 
with time. However, the discrete and the continuous 
equations are not equivalent when j ~ 0(1) because of 
the vanishing viscosity. 

Integrability is a remarkable feature of the Burgers 
equation: the Colc-Hopf transformation n = —vu x /u 
turns the nonlinear equation (fT]) into the ordinary dif- 
fusion equation u t = vu xx . Our discrete equation is 
not integrable. Interestingly, however, the discrete equa- 



tion 

= »i(»j- ( 31 ) 

which has a purely quadratic right-hand-side and mimics 
the nonlinear advection term in the Burgers equation, 
can be linearized. The transformation [2J] 

nj = (32) 
Qj+i 

reduces the nonlinear equation (|3ip into the linear equa- 
tion dqj/dt = qj-i- Using this transformation, equation 
(T5TT) can be solved exactly for the initial condition (|30l) . 
From the solution qj = 1 + t + jt 2 + ■ ■ ■ + 7jrjTrt J ~ 1 , 
the rarefaction wave discussed above follows immediately. 
However, the discrete equation (|3"Tj) does not satisfy con- 
tinuity. Constructing a discrete analog of the Burgers 
equation which satisfies the continuity condition, and in 
addition, can be transformed into a linear equation is an 
interesting challenge. 

In summary, we introduced a discrete nonlinear equa- 
tion that captures the key properties of the Burgers equa- 
tion. Intriguingly, the most natural discrete analog of 
the Burgers equation is essentially unique: matching the 
nonlinear advection term dictates that a viscosity that is 
equal to the density, v = n. The discrete equation satis- 
fies the continuity condition and consequently, the mass 
of an initially localized "bump" is conserved. The mass 
conservation ensures that the velocity of a traveling wave 
and the shape of a triangular wave match the continuous 
behavior, to leading order. We also find that the shape 
of a triangular wave is identical to the continuous case. 

In contrast with the classic Burgers equation, the vis- 
cosity equals the density. This discrepancy results in sub- 
tle departures from the continuous behavior: The shape 
of the traveling wave is asymmetric, and moreover, the 
front of the wave is extremely sharp with double expo- 
nential tail. Further, the front of the triangular wave 
includes a logarithmic-in-time correction to the leading 
behavior which is diffusive. 
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